Low Energy Excitations in Spin Glasses from Exact Ground States 



Matteo Palassini 

University of California, 3333 California Street, Suite 415, San Francisco, CA 94118 



m 
o 
o 

(N 

>>: 



C/3 



o 

> 

in 

(N 
O 



O 

o 



X 
S3 



Frauke Liers and Michael Juenger 
Institut fiir Informatik, Universitdt zu Koln, D50969 Cologne, Germany 



Physics Department, 



A. P. Young 
University of California, Santa Cruz CA 95064 
(Dated: February 1, 2008) 



We investigate the nature of the low-energy, large-scale excitations in the three-dimensional 
Edwards- Anderson Ising spin glass with Gaussian couplings and free boundary conditions, by study- 
ing the response of the ground state to a coupling-dependent perturbation introduced previously. 
The ground states are determined exactly for system sizes up to 12'' spins using a branch and cut 
algorithm. The data are consistent with a picture where the surface of the excitations is not space- 
filling, such as the droplet or the "TNT" picture, with only minimal corrections to scaling. When 
allowing for very large corrections to scaling, the data are also consistent with a picture with space- 
filling surfaces, such as replica symmetry breaking. The energy of the excitations scales with their 
size with a small exponent 9' , which is compatible with zero if we allow moderate corrections to 
scaling. We compare the results with data for periodic boundary conditions obtained with a genetic 
algorithm, and discuss the effects of different boundary conditions on corrections to scaling. Finally, 
we analyze the performance of our branch and cut algorithm, finding that it is correlated with the 
existence of large-scale, low-energy excitations. 

PACS numbers: PACS numbers: 75.50.Lk, 75.40.Mg, 05.50.-|-q 



I. INTRODUCTION 

There is still considerable debate about the nature of 
the spin glass state in finite dimensional spin glasses. 
Two principal theories have been investigated: the 
"droplet theory" proposed by Fisher and Huse^ (see also 
Refs. m, and the replica symmetry breaking picture 
of Parisi'*'^'^. In the droplet theory the lowest energy 
excitation of length scale I (a "droplet") has energy of 
order where is a positive exponent. Furthermore, 
the droplets have a surface with fractal dimension, ds, 
less than the space dimension d. 

Replica symmetry breaking (RSB) is well established 
in mean field theory, but it remains to be proven in finite 
dimensions. The precise nature of RSB in finite dimen- 
sions is not uniquely defined but it is generally agreed 
that a key feature of RSB is the existence of excitations 
whose energy, unlike that of droplets, remains of order 
unity even as their size tends to infinity. Furthermore, 
upon the creation of such a large scale, finite energy ex- 
citation, a finite fraction of the bonds change state (from 
satisfied to unsatisfied, or vice- versa) or, equivalently, the 
surface of these excitations is space filling, i.e. ds = d. 

Recently, Krzakala and Martin^ (KM) , and two of us® 
(PY), have argued, on the basis of numerical calculations 
at zero temperature, in favor of an intermediate scenario 
where there are large scale excitations whose energy does 
not increase with size, as in RSB, but which have a sur- 
face with ds < d. Following KM we shall denote this 
the "TNT" scenario. In the TNT scenario it is necessary 
to introduce two exponents which describe the growth 



of the energy of an excitation of scale L: (i) 9 (> 0) 
such that is the typical change in energy when the 
boundary conditions are changed, for example from pe- 
riodic to anti-periodic, in a system of size L, and (ii) 9' , 
which characterizes the energy of clusters excited within 
the system for afixed set of boundary conditions {9' was 
called 9g in Reffj). The TNT picture has been challenged 
(although in opposite senses) by Marinari and Parisi— and 
by MiddletoniSi. Subsequently, low temperature Monte 
Carlo simulations^ have found results consistent with 
the TNT scenario. The RSB, droplet, TNT and some 
other scenarios have been also studied by Newman and 
SteiniSii^. For some recent related work, see Refs-EHi 

The work of KM and PY determined the ground state 
with and without a certain perturbation (which was dif- 
ferent in the two cases), designed so that the ground 
state of the perturbed system is a large scale excitation 
of the original system. They used heuristic algorithms, 
i.e. algorithms which are not guaranteed to give the ex- 
act ground state, although both KM and PY argue that 
they do find the exact ground state in most cases. 

In this paper, we reconsider the problem of determin- 
ing 9' and ds, following the perturbation approach of 
PY, described in Section ^ but we apply an exact al- 
gorithm, known as "branch and cut" , so we are guar- 
anteed that the true ground state is reached every time. 
Exact optimization algorithms have been used before for 
spin glasses, see e.g. Refs. llTllTlTfll but, to our knowl- 
edge, their use in three-dimensions has been restricted to 
smaller sizes than studied here, and they were not used 
to investigate the real-space structure of the low-energy 
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excitations. 

Our implementation of the branch and cut technique 
can handle significantly larger sizes for free boundary 
conditions (be) than for periodic bo^li, so we use free 
be here. We consider a different (and enlarged) set of 
observables than PY, in the attempt to gain a fuller un- 
derstanding of what picture fits better the whole set of 
observables. We also perform a similar analysis of the 
data of PY, who used periodic be, in order to investi- 
gate the effects of different types of boundary conditions. 
The various pictures discussed refer to the large volume 
limit, while the sizes that can be currently reached are 
rather small. We will therefore pay particular attention 
to properly take into account corrections to scaling. In 
particular, we will try to determine what values of the 
parameters 9' and d — ds fit the data in the more "natu- 
ral" way, namely with the smallest corrections to scaling 
for the range of sizes considered. 

A summary of our results is as follows. We find that 
for periodic be, a simple scaling ansatz fits the results 
in a natural way, i.e. with negligible corrections to scal- 
ing and no adjustable parameters besides d — ds and 9' . 
This gives d ~ ds ^ 0.42 ± 0.03; 9' = -0.01 ± 0.03 (the 
meaning of the error bars will be explained later), which 
agrees with the results of PY, and is compatible with the 
TNT picture. We cannot rule out crossover to either the 
droplet or the RSB picture at length scales larger than 
our system sizes, but these scenarios, especially the lat- 
ter, would require larger corrections to scaling than the 
TNT picture. 

For free be, all forms of fitting require some corrections 
to scaling. The most natural scenario, in the sense ex- 
plained above, gives d-ds = 0.45±0.02; 9' = 0.18±0.03, 
with small corrections (of the order of 3%), which is 
compatible with the droplet picture. Allowing somewhat 
larger corrections (of order 10%), the data are also com- 
patible with 9' = 0, namely with the TNT picture. Fi- 
nally, if we allow for very large corrections, the data are 
also consistent with the RSB picture. 

In the second part of the paper, we analyze the perfor- 
mance of the branch and cut algorithm. We find that 
the number of elementary operations required to find 
the ground state increases exponentially with the size, 
as expected since computing a ground state of a three- 
dimensional spin glass system is an A/'P-hard problemSi. 
We also find, interestingly enough, that the CPU time is 
larger for samples in which there is an excited state close 
in energy to the ground state energy, yet different from 
the ground state in the orientation of a large number of 
spins. We are not aware of any previous quantitative 
measures of this trend, which we expect to be common 
to other algorithms as well. 

The rest of this paper is organized as follows. In Sec- 
tion ^1 we describe the method of perturbing the ground 
states to get information about low energy excitations, 
introduced by PY. Our results for the nature of the large 
scale, low energy excitations are given in Section IlIII A 
short description of the branch and cut algorithm used 



is given in Section IIVI and the performance of the al- 
gorithm is analyzed in Section [3 Our conclusions are 
summarized in Section IVll 

II. GROUND STATE PERTURBATION 
METHOD 

The Hamiltonian of the spin glass model is given by 
H = - ^JijSiSj, (1) 

where the sites i lie on a simple cubic lattice with N = 
spins in dimension d = 3, Si = ±1, and the Jy are 
nearest-neighbor interactions chosen from a Gaussian dis- 
tribution with zero mean and standard deviation unity. 
Free boundary conditions are applied in all directions. 

For a given set of bonds we determine the exact ground 
state using a branch and cut algorithm discussed in Sec- 
tion IIVI Let S*!*^-* be the ground state spin configuration. 
As in PY we then perturb the couplings Jij by an amount 

proportional to s'f^Sj'^^ in order to increase the energy 
of the ground state relative to the other states and there- 
fore possibly induce a change in the ground state. This 
perturbation, which depends upon a positive parameter 
e, is defined by 

where Nb — dL'^^^{L — 1) is the number of bonds in the 
Hamiltonian. We denote the unperturbed ground state 
energy by and the perturbed energy of the same 
state by Ei°\ The energy of the unperturbed ground 
state will thus increase exactly by an amount Ai?*^°) = 
s(o)-i;(o) =e. The energy of any other state, a say, will 
increase by the lesser amount AE^"^ = E^"^ — E^"^ = 

e qi^''°'\ where 9;^"'"'' is the "link overlap" between the 
states "0" and a, defined by 

ir' = ^^T.si''sl'\st's^;\ (3) 

in which the sum is over all the Nh nearest neighbor 
bonds. Note that the total energy of the states changes 
by an amount of order unity. 

As we apply the perturbation, the energy difference 
between a low energy excited state and the ground state 
decreases by the amount 

Ai;(o)-Ai;(") =e (1-9;^"-"^). (4) 

If there is at least one excited state such that E^"^ — 
< A£:(°) - A£;("), then one of these excited states 
will become the ground state of the perturbed Hamil- 
tonian. We denote the new ground state spin config- 
uration by and indicate by qi and g, with no in- 
dices, the link- and spin-overlap between the new and 
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old ground states 5) and , where q is defined as 

usual hyq = 1 / N J2 sl°'' S^°'^ ■ 

Due to the spin flip symmetry of the Hamiltonian l|T]l. 
the ground state is doubly degenerate, and therefore the 
distribution of q is symmetrioSS around q = 0. Hence, in 
the rest of the paper we will restrict ourselves to g > 
without loss of information. 

Consider the probability P(e, L) (with respect to the 
random couplings) that q is less than unity, i.e. that 
S'j-"^ and differ in a finite fraction of the spins. As 
discussed by PY, we assume that P(e, L) is dominated 

by those samples in which 5*1°^ and differ by flipping 
a single connected cluster of spins, with linear size L. 
Deviations from this assumption give rise to corrections 
to scaling, as pointed out by MiddletoniS, and will be 
analyzed in Section IIIII There are two energy scales in 
the problem: the typical energy above the ground state 
of such an excitation, which scales as (9' = 6 in the 
droplet picture), and the threshold energy of Eq. |@J, 
which scales as eL~('^~'^») since 1 — qi is proportional to 
the surface of the excitation, 1 — qi ^ j^-id-d^) ^ Hence, 
the dimensionless probability P(e, L) is a function of the 
ratio of these two energy scales: 

Pie,L)^gieL-n, (5) 

where g{x) is a scaling function and 

^1 = 9' + d-ds. (6) 

From this we obtain scaling relations for various observ- 
ables. For example, since 1 — q ~ 0(1) and 1 — ~ 
L-(d^d,)^ we obtain^: 

(1-9) = F,{eL-^^) (7) 
{l-qi) = L-'^'-'^^^F.^ieL-n, (8) 

where (• ■ ■) is the average with respect to the random 
couplings. By measuring (1 — q) and (1 — qi) we can then 
determine d—dg and 9' , the two exponents discriminating 
the various pictures of the spin glass phase discussed in 
Section m 

For small e, we expect the probability that the ground 
state changes to be proportional to e (for fixed L) , which 
implies g{x) ~ a: for a; ^ 0. Hence Fq{x) and Fq^{x) 
also vary linearly for small x, and the asymptotic scaling 
behavior for L ^ e^/^ is 

(1-g) ^ eL-^\ (9) 

(l-qi) ^ eL-^\ (10) 

where 

^ii = 9' + 2{d - ds) . (11) 

In the RSB case, d — dg = 9' = 0, and therefore /i = 
Hi = 0. The scaling relations in Eqs. (|7||Hll reduce in this 
case to 



and the asymptotic scaling behavior to 

(l-q)^e, (l-<z,>^e (RSB). (13) 

We see that both scaling and asymptotic scaling are in 
a sense trivial in RSB since the L dependence disappears. 
Nevertheless, we will still use the term scaling. 

It is also convenient to analyze just those samples in 
which the unperturbed and perturbed ground states are 
very different, i.e. where q < gmax, a threshold value. 
Denoting such restricted averages by (•••)£, we have 

{l-qi)c = L-(''-'^^F^^{eL-n. (14) 

This is of the same form as in Eq. ijSJ), but, for sufh- 
ciently small Qmax, the behavior of the scaling functions 
Fq^ (x) and F^^ (x) at small argument will be different for 
the following reason. If we average over all samples we 
need to include the probability P{e, L) that the pertur- 
bation generates an excitation with q < 1. This is pro- 
portional to eL"'' for eL~'^ ^ 1, which is the reason why 
Fq{x) ^ X for small x. However, this factor is automat- 
ically taken into account in the selection of the samples 
in the restricted average in Eq. (|14|l . and so should not 
be included again when performing the average. As a 
result, Fq^ (x) tends to a constant for a; — > 0, therefore 
the asymptotic scaling is 

(l-qi),^ L-<^''-'^\ (15) 

In particular, in RSB this becomes 

(1 - (7,)c const. (RSB). (16) 

Note that in both cases the asymptotic scaling is inde- 
pendent of e. 

When analyzing the numerical data, we must be aware 
that there are corrections to both (simple) scaling and 
asymptotic scaling that occur when L is not large enough. 
Corrections to simple scaling take the form of additive 
corrections to relations such as Eqs. ||SJ), 0, ©, and 
(|14ll , whose amplitude is characterized by a correction-to- 
scaling exponent w. For example, including the leading 
correction, Eq. (|14|l becomes 

(1 - qi)c = [piieL-'^) + l-Gq,{eL-n] • (17) 

For eL~'^ 0, this gives the correction to asymptotic 
scaling corresponding to Eq. H15|l 

(^-^i)c = jj^[^+^)- (18) 
For the RSB case, this goes over to 

(l-qi), = a+A, (19) 



{l-q)^Fq{e), {1 - qi) ^ Fq,{e) (RSB), (12) 



rather than Eq. (|16|l . 
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' 4 




e/r = 1 


e/r = 2 


e/r = 4 


4 


5nnnn 


5nnnn 


50000 


50000 


50000 


6 


20000 


20000 


20000 


20000 


20000 


8 


15000 


13467 


13467 


6000 


6000 


10 


10000 


7440 


6000 


4918 


4000 


12 


5670 




4202 







TABLE I: Number of independent realizations of the disorder 
(samples) used in the computations. 



Even when these corrections to (simple) scaling are 
negligible and the scaling /orm, such as Eq. (|14|l . is valid, 
the argument of the scaling function may not be suf- 
ficiently small for asymptotic scaling to hold. In this 
regime, when fitting the data to asymptotic scaling we 
have to consider further corrections to (asymptotic) scal- 
ing, whose form is obtained by expanding the scaling 
function in its argument. For example, the leading cor- 
rection to Eq. (|15|l . coming from expanding the F^^ in 
Eq. 114|) to second order, will be 



(1 - qi)c = 



1 



(20) 



which goes over to (1 — qi)c = a + 5e in RSB. In general, 
both types of corrections need to be borne in mind when 
fitting the data. 



III. RESULTS 

We applied the perturbation method described in the 
previous Section to systems of size L = 4,6,8,10, and 
12. For each size, we considered five values of the pertur- 
bation strength e, namely e/r = i, i, 1,2, and 4, where 
T — \/& is the mean field transition temperature, except 
for i = 12 for which only e/r — ^ and 1 were considered. 
We choose this value of r so we can compare our results 
with the results of FY for periodic be. In order to dis- 
criminate between the different pictures, it is important 
to have high statistics. Table reports the number of 
samples computed for each size. Note that the number 
of samples necessary to achieve a given statistical error 
increases as e decreases, since the fraction of samples in 
which ^ 3'^°^ decreases. 



A. Spin and link overlap 

1. Qualitative analysis 

We start with a qualitative analysis of the results. In 
Fig. n we show scatter plots in the (g, qi) plane for L — 
4, 10 and £/t — i, 1, where each point represents one of 
2000 randomly generated samples. Clearly, the link- and 
spin-overlap are strongly correlated. We note that, as e 
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FIG. 1: Scatter plots in the {q,qi) plane for different values 
of the size L and perturbation strength e. 



decreases, there are less points with small g, and that 
as L increases the data shift towards larger values of qi . 
Similar plots^^ for periodic be show that qi is significantly 
lower than for free be. While q has a large variance (the 
points are distributed along the whole q axis), the link- 
overlap has a much smaller variance, which decreases as 
L increases, suggesting that in the thermodynamic limit 
qi either tends to one or becomes a well defined function 
of q. 

To quantify this, in Fig. [21 we show the standard devi- 
ation of qi 



(21) 



restricted over samples with q < qmax, as a function of 
the system size for e/r = 1. We take Qmax — 0.2, since 
we are interested in the region of small q, which corre- 
sponds to large-scale excitations. A power law a = aL~^ 
fits well the data with S = 0.52 ± 0.03 (x^ = 1-80, the 
best fit is shown in Fig. O. Here and in the following, 
unless stated otherwise, the error bars on the fit param- 
eters are purely statistical in relation to the fitting form 
considered^fi-. Restricting the average in Eq. H21|l to dif- 
ferent intervals of q gives results also compatible with a 
power law. A vanishing a in the thermodynamic limit 
is consistent with RSB, which predicts that qi is a (non- 
trivial) function of q. It is also trivially consistent with 
the droplet model or the TNT picture, where qi ~ I for 
all q. 

We also measured how the standard deviation of q 
varies with L, finding that it varies between 0.28 and 
0.32. It can be fitted both to a constant (as expected in 
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q|< = 0.2 



e = V 6 

4 5 



6 7 
L 



8 9 10 




FIG. 2: A plot of the standard deviation of the link-overlap 
= {{if}': ~ {li}^V^^ J where the average is restricted to sam- 
ples such that q < 0.2. The line represents a power law fit 
with exponent S — 0.52. 



FIG. 3: Average of the link-overlap restricted to intervals of q 
of width 0.1, as a function of q for different sizes L (from bot- 
tom to top, L = 4,6, 8, 10, 12). The continuous lines represent 
quadratic fits including values of q up to where the lines end. 
The dashed lines are a guide to the eye. 



RSB) or to a power law with a small exponent around 
0.1. However the error bars are very large hence the fits 
are not very informative. 

Under the RSB hypothesis, it is interesting to study 
the functional relationship between q and qi. In Fig. |3 
we show the average value of qi, restricted to intervals 
q G [(Zmin, 9max], as a function of the mean value of q 
in each interval^i for e/r = 1. For fixed L, a quadratic 
form qi = a{L)+ P{L) q^, motivated by the infinite-range 
Sherrington-Kirkpatrick model where qi = q^, fits well 
the data for q less than some cut-off value which increases 
with L (see Fig.l^J. The quadratic fit works well also for 
other values of e, and a{L) and /3{L) show a modest 
variation with e. We also tried global fits including data 
for all values of q and L, obtaining similar results. 

Extrapolating a{L) and P{L) to L ^ oo with fits of 
the form a{L) = a + b/L", (3{L) =13 + b' /L"' , we obtain 

qi = (0.77 ± 0.02) + (0.27 ± 0.03) q^ , (22) 

where again the errors are purely statistical for the func- 
tional form considered. This nontrivial relation between 
q and qi in the thermodynamic limit is consistent with 
RSB, while in the droplet or TNT pictures the data in 
Figl^l would shift to = 1 in this limit. 

The power law form 1 - a{L) = b/L", (3{L) = b' / L"' , 
which implies = 1 in the large volume limit, fits poorly 
the data if we include all sizes. However, if we exclude 
the L = 4 data, the quality of the fit becomes as good 
as that of the RSB fit just discussed. Hence, allowing for 



small corrections to scaling, the droplet or TNT scenario 
are consistent with the data. 

This already shows that care must be taken to prop- 
erly consider corrections to scaling when comparing the 
merits of the fits to various pictures. In the following, 
we will investigate in detail the validity of the various 
pictures by considering several observables and explicitly 
discussing the corrections to scaling for each picture. 

2. Determination of d — da 

We start with the determination of d — from various 
observables. We will show that for all observables, a wide 
range of values oi d — dg fits well the data when allowing 
corrections to scaling, but that for all observables the 
smallest corrections are attained for a value of d — c?s 
around 0.44 as in PY. 

The main part of Fig. ^ plots (1 — qi)^ as a function 
of L, for various values of e and gmax = 0.4 (^max = 0.2 
gives essentially the same results). First, note that, inde- 
pendently of what picture holds in the L — > oo limit, the 
data deviate significantly from asymptotic scaling, see 
Eq. (|15|) . in which the various e values should collapse 
on a single curve. Second, the data have a noticeable 
positive (upward) curvature for all values of e. In Sec- 
tion llll Cl we will show that in the case of periodic be the 
data have a much smaller dependence on e and a much 
smaller curvature (see Fig. I12|l . 
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In order to determine how the various pictures fit the 
data of Fig. ^ we start by considering, following Ref|3, 
the following three functional forms: 



free b.c. 



Form I 
Form II 
Form III 



{l-qi)c 
{l-qi)c 
{l-qi)c 



a + b/L" 

a + b/L + c/L^ (23) 
b/L" 



Form I corresponds to the RSB prediction including 
the leading correction to scaling, see Eq. fH^ . with c = u. 
Form II is a different parameterization of the corrections 
to scaling. Form III corresponds to the asymptotic be- 
havior of both the TNT and droplet pictures without cor- 
rections to scaling, see Ec^. (|15|l . with c = d — dg- 

The results of these fits (performed by minimiza- 
tion) are reported in Table m From the Table we see 
that Forms I and II, appropriate to RSB, fit well the 
data with a low and a > outside the error bars. The 
best fits with Form I are shown by the dashed lines in 
Fig. 2] The variation of a between Forms I and II pro- 
vides a measure of the systematic error associated with 
the unknown corrections to scaling. Within this error, 
a is independent of e, as predicted by RSB. Therefore, 
the data for (1 — qi)c are compatible with RSB, and our 
central estimate under the RSB hypothesis is 



lim (1 - qi)c. = 0.20 ± 0.02 



(RSB) , (24) 



where the error takes into account also the uncertainty in 
the form of the corrections to scaling, assuming that the 
corrections considered in either Form I or II describe well 
the data in the whole range of sizes considered. Marinari 
and Parisi^ fitted Form I (resp. II) to their data for peri- 
odic boundary conditions, L < 14, and e/r = 4, and ob- 
tained a = 0.24 (resp. a = 0.30), from which we estimate 
a central value a = 0.27 ± 0.05. This is just in agreement 
with our estimate above for free be, suggesting that the 
infinite volume limit of {l^qi)c, if nonzero, may be inde- 
pendent of the boundary conditions, although we do not 
have an argument why this should be the case. 

The power law Form, III, appropriate to the droplet 
model or the TNT scenario, does not fit well the data if 
we include all the sizes, but if we exclude L = 4, it does 
fit well for e/r < 2, and the fit parameters b and c are 
almost independent of e. The quality of the fit of Form 
III is worse than that of Forms I and II, but still accept- 
able. The main point we want to stress, however, is that 
the worse fit of Form III alone does not necessarily favor 
the RSB picture, since Form III does not include correc- 
tions to scaling, while Forms I and II do. In other words. 
Forms I and II are rather "forgiving" with the RSB pic- 
ture, allowing corrections of magnitude 100% - 200% of 
the predicted L-independent asymptote, while Form III 
demands that the power-law scenario fits with corrections 
smaller than the (very small) statistical errors. By look- 
ing at Figure 01 it is clear that the data are closer to a 
power law than to an L-independent behaviour. 

Therefore, in order to try a comparison that puts the 
various pictures on an equal footing, we performed fits 




FIG. 4: Logarithmic plot of the average (1 — qi}c, restricted 
to samples with |q| < 0.4, as a function of the system size 
L. Only three values of e are displayed for clarity. The lower 
continuous line is the best fit with a power-law. Form III in 
Eq. 1231 for e/r = j, where r — y/6, and the L = 4 data have 
been excluded from the fit. The dashed lines are the best fits 
with Form I in Eq. I|28|l . The inset shows a scaling plot of the 
data in the main figure, excluding the L = 4 data, according 
to Eq. (1141 . Here and in the following figures, note that the 
data for various e are correlated, since the samples used for 
large e were also used for small e. 



with the following more general functional form: 



Form IV: 



-d-d. 



{l-qi)c = a + b/L' 



(25) 



where we fix d — ds and minimize the with respect to 
a, b, c, repeating the procedure for different values of d — 
ds- For d — ds = 0, Form IV reduces to Form I, while for 
d — ds > 0, it corresponds to Form III plus a correction- 
to-scaling term, with correction-to-scaling exponent lu = 
c. We find that, as we might have expected from the 
previous discussion. Form IV fits well the data for a wide 
range of values of d — dg. For example, for e/r = i, a 
value of d — ds between and 0.45 gives a goodness-of-fit 
parameter Q > 0.43, which is entirely acceptable. 

This shows that, when allowing for corrections to scal- 
ing for all pictures, the droplet or TNT pictures are as 
good as RSB as far as the statistical quality of the fits 
is concerned. However, within the interval of acceptable 
values of d — ds, clearly the larger is d — ds the smaller 
are the corrections to asymptotic scaling. For example, 
for e/r = j and d — ds = 0.42, the correction term b/L'^ 
in Form IV amounts to only 0.1% of the total for L = 12, 
while for d — ds = it amounts to 43%. Hence a large 
value of d — ds may be regarded as more "natural" in this 
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Form 


./r 




Q 


a 


b 


c 




0.25 


0.014 


0.99 


0.171(2) 


1.068(5) 


0.890(8) 




0.5 


0.015 


0.90 


0.185(2) 


1.14(1) 


0.96(1) 


I 


1 


3.04 


0.22 


0.201(14) 


1.26(7) 


1.04(7) 




2 


0.39 


0.52 


0.206(6) 


1.42(3) 


1.08(3) 




4 


0.27 


0.60 


0.215(4) 


1.73(3) 


1.14(2) 




0.25 


0.037 


0.98 


0.182(1) 


1.29(2) 


-0.40(5) 




0.5 


0.010 


0.92 


0.189(1) 


1.23(1) 


-0.14(3) 


II 


1 


3.01 


0.22 


0.198(8) 


1.16(10) 


0.16(27) 




2 


0.47 


0.49 


0.199(4) 


1.21(5) 


0.31(14) 




4 


0.42 


0.52 


0.202(3) 


1.31(4) 


0.64(11) 




0.25 


1.40 


0.49 




0.85(2) 


0.44(1) 




0.5 


1.63 


0.20 




0.87(3) 


0.45(2) 


III 


1 


9.81 


0.007 




0.88(3) 


0.44(2) 




2 


6.75 


0.009 




0.95(4) 


0.47(2) 




4 


11.4 


0.0007 




1.09(5) 


0.51(2) 




0.25 


0.55 


0.76 


0.808(4) 


7(8) 


3.6(8) 




0.5 


0.65 


0.42 


0.811(6) 


6(7) 


3.2(8) 


IV 


1 


5.91 


0.05 


0.828(7) 


18(37) 


4.0(1.5) 




2 


2.80 


0.09 


0.844(9) 


5(5) 


2.9(7) 




4 


3.01 


0.08 


0.87(1) 


3.8(1.6) 


2.3(3) 



TABLE II: Fits to (1 — qi)c with ijmin = and (jmax = 0.4. 
The four groups of data refer, from top to bottom, to the 
three fitting functions I, II, and III in Eq. 1)23^ . and Form IV 
in Eq. (O with d- ds = 0.44. For Form III, data for L = 4 
was not included in the fit. The errors are the standard errors 
of a nonlinear fitting routine—, and Q is the goodness-of-fit 
parameter. 

range of sizes. 

If we impose that the correction to scaling for L > 8 
is less than an (arbitrary) limit of 3%, we obtain the 
estimate 

d- 4 = 0.44 ± 0.03 (26) 

where the error is purely statistical within this assump- 
tion. In Table II we show the fits obtained with Form 
IV imposing this value. This agrees with the estimate 
d- dg = 0.42 ± 0.02 of PY for periodic be (see also Sec- 
tion I^^C] of this paper), which is reassuring since d — dg 
should not depend on the boundary conditions. 

Note that for d — dg =0, corresponding to Form I, 
corrections within 3% from the asymptotic limit would 
only be attained for a size L ~ 200. We also note that, 
as we discussed in Section ^ even in the regime where 
corrections to scaling are negligible, asymptotic scaling 
sets in only for L ^ e^/^. This explains why, if d — 
dg ~ 0.44, the quality of the power-law fit in Table II 
becomes progressively worse as e increases. In particular, 
the deviation from asymptotic scaling is very significant 
for e/r = 4, and hence from the data of e/r = 4 alone 
one should not necessarily concludeSi that an asymptotic 
power-law behavior is ruled out. This is seen also in the 



inset of Fig. 01 which shows that, if we exclude L = 4, 
the data are compatible with the scaling relation Eq. ||SJ), 
where the exponent is independently determined below. 

PY determined d—dg from the ratio R = (1 — — g) 
which has the same scaling behavior as the quantity 
(1 - qi)c used here, namely R = L-^'^-'^''^ FR{e/ L^^), with 
Fr{x) ^ const, as a; ^ 0. Middleton^" observed that, in 
two dimensions, small droplets introduce significant cor- 
rections to scaling, and suggested that this may be the 
case also in three dimensions, possibly invalidating the 
conclusions of PY. The quantity {I — qi)c is less sensitive 
to these corrections since, by restricting to small q, small 
droplets should give a smaller contribution, because we 
have eliminated the part at large q where one can have 
only small droplets. Hence, to investigate these correc- 
tions, we fitted our data for R with Forms I-IV above 
(with R replacing (1 — g;)c). The results we find are very 
similar to those for (1 — qi)c: Forms I and II fit well the 
data with a low x^, giving a = 0.27 ±0.03 independent of 
e within the error bars. A power law fits well the data if 
we exclude L = A, with an exponent d — dg = 0.43 ± 0.03 
nearly independent of e and in agreement with Eq. H2t)|l . 
The residual dependence on e is well accounted for by a 
scaling plot similar to the inset in Fig.^J Form IV also fits 
well the data for a wide range of values of d — dg. Again, 
a power law is more natural in the sense that corrections 
to scaling are smaller, and the smallest corrections are 
obtained for d — dg around 0.43 as for (1 — qi)c- We in- 
terpret the fact that the two quantities give the same 
value oi d — dg as evidence that corrections due to small 
droplets are indeed not important in three dimensions in 
this range of sizes. In Section llll Cl we will show that this 
is also the case for periodic be. 

To summarize this part, the data for both R and (1 — 
qi)c are compatible with a wide range of values of d — dg 
between zero (corresponding to RSB) and ~ 0.44, but a 
value at the higher end of this range describes the data 
in a more natural way, in the sense that the corrections 
to scaling are smaller. 

3. Determination of 9' 

Next, we turn to the exponent defined in Eq. lO, 
from which we will extract the exponent 9' which is the 
other exponent, with d—dg, characterizing the spin glass 
phase. To this end we consider the ratio 



which follows the scaling law 

B^FBie/Lf^). (28) 

The factor L'^~'^' does not appear here since it cancels 
between numerator and denumerator of Eq. H27|l . thus 
allowing us to determine ^ independently oi d — dg. Fol- 
lowing the arguments in Section^ we expect Fb{x) ^ x 
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for small x since both L'^''^^ (l-qi) and L^Cd-d.) ^(i_^^)2^ 
vary as e/L^ for e/L'' 0, hence the asymptotic scaling 
of B is B - e/L^. 

To determine fi, wc fit the scaling law Eq. (|28|l to our 
data assuming a polynomial form of order n for Fb{x), 
namely Fb{x) = J2i=o n^^^^ ^ with co set to zero in order 
to satisfy the asymptotic scaling Fb {x) 0. We 

repeat the fit in an interval of values for /i, and determine 
the value of /i which gives the best fit, varying n until 
the of the best fit becomes approximately constant. 
In this way we obtain 



free b.c. 



= 0.63 ±0.03, 



(29) 



where the error is purely statistical, under the assump- 
tion that the corrections to scaling are smaller than the 
statistical errors of the data. As shown in Fig. |S1 scal- 
ing is quite satisfactory, with all the data collapsing on 
one curve, although the data for different e overlap only 
slightly. The best fit for n = 6 is displayed by the con- 
tinuous line. We emphasize that this scaling plot is ob- 
tained with only one adjustable parameter, fi. Note that 
in the asymptotic scaling regime the data should follow 
a straight line (power-law), while the data in the figure 
show a pronounced curvature. Significant corrections to 
asymptotic scaling must be expected for large e, since B 
must satisfy the inequality B < 1. The dashed line in 
Fig.Elrepresents the linear term of Fb{x), corresponding 
to asymptotic scaling, and the deviation from it gives a 
measure of the corrections to asymptotic scaling. The 
e/r = i data are quite close to asymptotic scaling, while 
the data for large e deviate significantly from it. An- 
other manifestation of these corrections is that, if we fit 
the data with a power-law B = b/L'^, the effective ex- 
ponent jl varies strongly with e, converging towards 0.63 
for e ^ 0. 

In RSB, i? ~ e as L ^ oo since ^ — 0. To test the 
RSB prediction, we performed fits of B using Form I 
and Form II in Eq. (|23|l . (where (1 — g;)c is replaced by 
B). Form I gives unphysical (negative) values of a, while 
Form II gives an acceptable fit with a positive a roughly 
proportional to e. Therefore, the data for B cannot rule 
out RSB. Note that, if RSB holds asymptotically, the 
data in Fig|Sl would deviate from the scaling curve for 
larger L, saturating to a constant value for small values 
of e/L^. The good data collapse we observed, therefore, 
would be entirely accidental. 

We believe that the observed data collapse is a good 
indication towards the validity of a scaling scenario with 
large fj,. Certainly this scenario is more natural since 
it fits the data with (almost) no corrections to (simple) 
scaling, while the corrections for RSB are very large as 
apparent from FigEl A similar conclusion was reached 
in the determination of d — ds- 

As a further test, we can obtain a second estimate 
of fi from the quantity {1 — qi) {q unrestricted), whose 
scaling and asymptotic scaling are given in Eqs. © and 
(|10|l . We find results similar to those for B: A power 
law fit (1 — qi) — b/Li^' (Form III) gives acceptable fits 
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FIG. 5: Scaling plot of the ratio B = (1 - gi)V((l ^ nf) 
according to Eq. 1281 1. The continuous line is a polynomial fit 
of order n = 6, which gives x^/d.o.f = 1.09, and a goodness- 
of-fit parameter Q — 0.35. The dashed line is the linear term 
of the polynomial fit, corresponding to the asymptotic scaling 
for L oo. 



for e/r < i and for all values of e if we exclude L = 4. 
As for B, the effective exponent /!/ changes significantly 
with e, due to corrections to asymptotic scaling, and by 
extrapolating it to e = we obtain /i; = 1.15±0.08. This 
gives /i = /i/ — (d — ds) = 0.60 ± 0.09 which agrees with 
the estimate /i — 0.63 ± 0.03 obtained from B. We also 
verified that, as for B, the data collapse reasonably well 
on one curve for /i = 0.64 ± 0.05 according to Eq. ((SJ, 
although the quality of the scaling is somewhat worse 
than that of Fig[31 To check the RSB prediction, we 
fitted the data to Forms I and II (where now (1 — qi)c 
is replaced by (1 — qi)), finding that they both fit well 
the data, with a roughly proportional to e as expected 
in RSB, although, for small e, a is also compatible with 
zero. Therefore, as for B, the data are also consistent 
with RSB, but this scenario requires large corrections to 
scaling, while the hypothesis fi — 0.63 fits the data with 
almost no corrections to (simple) scaling. 

In the analysis so far, we have determined the expo- 
nents fi and d — ds using just the link-overlap qi. By con- 
trast, PY determined /it (for periodic be) from the scahng 
of the spin-overlap q. An advantage of qi is that its vari- 
ance is much lower, as shown in Fig. ^ In any event, we 
have verified that the scaling relation Eq. Q fits well the 
data for g, giving ^ = 0.65 ± 0.02, in agreement with the 
estimates from B and {1 — qi). 

Summarizing this part, we find that the data for all the 
quantities considered, namely B, (l — qi), and (l^q), are 
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consistent with the RSB prediction that /i = asymptot- 
ically, but large corrections to scaling are required in the 
fit, similarly to what we observed in the determination of 
d—dg. The data are also fitted very well by a scaling sce- 
nario with /i ~ 0.63, with almost negligible corrections 
to scaling (but with sizeable corrections to asymptotic 
scaling, which instead were small for the observables con- 
sidered for d — dg). Under the "natural" assumption of 
small corrections to scaling, from the estimates of ^ and 
d~ ds in Eqs. and (j^ . we obtain 



Ai- (d-d,) = 0.19 ±0.06, 



(30) 



where, again, the error is purely statistical subject to 
the condition of having small (less than 3%) corrections 
to scaling. This result agrees with the droplet theory 
which predicts that 9' — 9 > 0, and is compatible with 
the value of 6* ~ 0.2 characterizing the energy of domain 
walls induced by a change in boundary conditionaiSi. By 
contrast, for periodic be and under the same assumption 
of small corrections to scaling, 9' is compatible with zero 
(see PY and Section UlI Cll . In Section UlIDI we will an- 
alyze the origin of this discrepancy, and show that by 
allowing small (of order 10%) corrections to scaling the 
free be data can be reconciled with 9' ~ 0. 



B. Box overlaps 

So far we have analyzed the link- and spin-overlap 
which are computed on the whole system (bulk). We now 
turn to a different observable, the box-overlap defined as 



IB 



Td 1 



(0) ^(0) 



(31) 



where the sum runs over the sites contained in a central 
cubic box of fixed size Lb = "2. In the following we will 
only consider the absolute value {qsly which we still call 
qb for simplicity. When a large-scale cluster of spins is 
flipped, for large L the probability that its surface goes 
across the central box is proportional to the ratio of its 
surface area, ^ to the volume, L'^. Therefore 1 — 
qB ~ L^^'^^'^^'i from which we obtain the scaling laws 



(32) 
(33) 



where, as for the corresponding scaling functions for qi, 
Fq^{x) r^x imd F^^{x) - 
asymptotic scaling for L 



const, for small x. Hence the 
-!■ oo is 



(l-gs) - eL-^' 



(34) 
(35) 



In RSB, this reduces to (1— gs) ~ e and (1— g_B)c ~ const. 
An advantage of qs over qi is that the former, being 
measured away from the boundaries, should have smaller 
corrections to scaling and be less sensitive to boundary 




4 5 



6 7 8 910 
L 



FIG. 6: Logarithmic plot of the average box-overlap, re- 
stricted to samples such that q < 0.4. We show the data 
for just two values of e for clarity. The data for other values 
of e are superimposed. The lower continuous line is a power- 
law fit for e/r = 4. The dashed line is the fit with Form II in 
Eq. 123L with qb replacing qi. The slope gives the exponent 
d — ds- 



conditions. Indeed, Monte Carlo simulations^^-^^ show 
that qB has rather small corrections to scaling. This is 
likely to be particularly important for the free boundary 
conditions used here. 

Fig. |B| shows the restricted average (1 — '?b)c, with 
qmax = 0.4, as a function of L for two values of e. The 
data are clearly decreasing with L, are essentially inde- 
pendent of e, as expected from Eq. (|35|l . and are close 
to a straight line on the logarithmic plot. This indi- 
cates that the power law fit. Form III, appropriate to the 
droplet and TNT scenarios, should work well and indeed 
it does, even for the largest value of e (we note however 
that the statistical errors are larger than for the link- 
overlap, hence the fits are less sensitive to corrections to 
scaling). The exponent is almost independent of e, vary- 
ing between 0.48 and 0.52, and from this we obtain the 
estimate 



d-d, = 0.48 ±0.03 



(36) 



which is in agreement with the estimates d~ ds — 0.44 ± 
0.03 from (1 - qi)^ and d-ds= 0.43 ± 0.03 from R. 

Forms I and II (with qB replacing qi) also fit well the 
data, with a between 0.14 and 0.36 (with no discernible 
trend with e). Hence the data are also compatible with 
RSB, and under the RSB hypothesis, we estimate 



lim (l-gi3)c = 0.25 ±0.10 (RSB). 



(37) 
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FIG. 7: Logarithmic plot of the average box-overlap, mul- 
tiplied by r/e in order to highlight the deviation from the 
asymptotic behavior of Eq. 13411 in which the data for various 
e should collapse on a single curve. The continuous lines rep- 
resent fits with the power-law Form III excluding L = 4. The 
dashed lines represent fits with Form I in Eq. (I23II . 



FIG. 8: Scaling plot of the box-overlap according to Eq. (1321 . 
The continuous line is a polynomial fit of order n = 6, which 
gives x^/d.o.f = 0.63, and a goodness-of-fit parameter Q = 
0.85. The dashed line is the linear term of the polynomial fit, 
corresponding to the asymptotic behavior for L ^ oo. 



As usual, we note that the RSB scenario requires rather 
large corrections to scaling, while the power law fits the 
data with no corrections. 

Fig. [3 shows the unrestricted average (1 — qb) multi- 
phed by r/e, which asymptotically should be indepen- 
dent of e. The data show a small curvature and a signif- 
icant e dependence, indicating that for this quantity we 
are not yet in the asymptotic scaling regime (similarly to 
what we observed for the quantity B). The data are fit- 
ted well by a power law, with an exponent that changes 
with e and tends towards /x ~ 0.63 for e — > 0. Fits using 
Forms I and II give a compatible with zero. We also de- 
termined /i from the scaling relation Eq. H32|l . by fixing 
d — ds — 0.44 and using the same fitting procedure as for 
B (which assumes no corrections to scaling), finding 

^ = 0.62 ±0.04 (38) 

which agrees with the various estimates of fi obtained 
from B, {1 ~ qi), and (1 — q). Fig. |H1 shows the corre- 
sponding scaling plot, in which the data collapse is fairly 
good. 

To conclude this subsection, the data for box overlaps 
can be fitted with smaller corrections to scaling than the 
data for the bulk link- and spin-overlap. A fit to the 
generic scaling picture, with no corrections to scaling, 
gives results for the exponents d — dg and fj, in agreement 
with those from the bulk quantities analyzed in the previ- 
ous subsections. However, as with the bulk observables, 



assuming large corrections to scaling, the data can also 
be fitted to the RSB picture. 

C. Comparison with periodic boundary conditions 

In order to assess the effect of different boundary con- 
ditions, we have repeated part of the analysis above (with 
the exclusion of box-overlaps) for the data of FY (Ref'^ 
for periodic boundary conditions and L < 8. The ground 
states were obtained using a hybrid genetic algorithm as 
described in FY. This does not guarantee to find the true 
ground state, but the systematic errors due to occasion- 
ally missing it arc smaller than the statistical errorA 

Figs. |51 and ^] show the equivalent for periodic be of 
Figs. 01 and for free be. The data in Fig.lHlshows much 
less curvature and also a smaller dependence on e than 
the corresponding data for free be in Fig. 01 indicating 
that corrections to scaling are smaller for periodic be. 
Table UTTI reports the best fits using the three functional 
forms of Eq. (|23|l . Form I fits well the data, but a varies 
significantly with e, and for small e it is compatible with 
zero. Form II also fits well, with a independent of e within 
the statistical errors. From this fit we estimate 

lim (1 - = 0.28 ±0.03 (RSB) (39) 

L — ^oo 

(see the comment after Eq. (|24ll as to the meaning of the 
error bar) which agrees with the estimate 0.24 of Mari- 
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FIG. 9: Same as Fig.|l]but for periodic boundary conditions, 
using the data of PY (RefiQ). 



FIG. 10: Same as Fig.|^but for periodic boundary conditions, 
using the data of PY (Ref|g). 



nari and ParisiS-, and is just consistent with our estimate 
0.20 ± 0.02 for free be. 

The power-law fit with no corrections, Form III, fits 
well the data for the two smallest values of e and, if we 
exclude L = 3, for all but the largest value of e. The 
exponent c = d — ds is nearly independent of e and gives 

d-ds = 0.43 ± 0.02 (periodic be) . (40) 

This result agrees with the estimate d — dg = 0.42 ± 0.02 
of PY obtained from the ratio R defined above, confirm- 
ing that corrections due to small droplets should not be 
important in three dimensions, and with our estimate 
d — ds = 0.44 ± 0.03 for free be, indicating that d — dg 
does not depend on boundary conditions. 

We also performed fits with Form IV which includes 
corrections to scaling. As for free be, a wide range of 
values of d — ds from zero to around 0.44 give a good fit, 
with the largest values giving the smallest corrections to 
scaling. The results of the fit for d — ds = 0.43 are shown 
in Table IIIII For the two smaller values of e, the fits 
are difficult because corrections to scaling are very small, 
hence they are not shown. 

We determined the exponent /i from the ratio B using 
Eq. (|28|l and the fitting procedure described for free be, 
obtaining 

^ = 0.42 ± 0.03 (periodic be) (41) 

where, as for the estimate of d — dg above, the errors are 
purely statistical with the assumption that corrections to 
scaling are smaller than the statistical errors of the data. 
Scaling is rather satisfactory as shown in Fig. ^| This 



Form 


e/r 




Q 


a 


b 


c 




0.25 


0.003 


0.99 


-0.076(7) 


1.256(6) 


0.384(4) 




0.5 


0.92 


0.62 


0.05(4) 


1.16(2) 


0.47(3) 


I 


1 


1.58 


0.45 


0.10(3) 


1.16(2) 


0.52(3) 




2 


2.20 


0.33 


0.12(3) 


1.18(2) 


0.54(4) 




4 


1.58 


0.45 


0.20(2) 


1.270(4) 


0.68(2) 




0.25 


0.33 


0.56 


0.279(7) 


1.90(5) 


-1.5(1) 




0.5 


5.22 


0.073 


0.28(1) 


1.90(9) 


-1.5(2) 


II 


1 


0.69 


0.71 


0.280(4) 


1.89(3) 


-1.40(6) 




2 


0.04 


0.98 


0.283(1) 


1.90(1) 


-1.33(2) 




4 


0.36 


0.83 


0.291(2) 


1.86(2) 


-1.01(4) 




0.25 


0.02 


0.887 




1.204(2) 


0.433(1) 




0.5 


0.35 


0.838 




1.193(3) 


0.427(2) 


III 


1 


5.14 


0.076 




1.21(2) 


0.434(6) 




2 


7.82 


0.020 




1.24(2) 


0.440(8) 




4 


25.7 


2 lO"'^ 




1.31(2) 


0.46(1) 




1 


3.59 


0.16 


1.205(4) 


0.3(1.5) 


3(4) 


IV 


2 


4.69 


0.09 


1.214(7) 


0.2(5) 


2(2) 




4 


8.38 


0.01 


1.231(8) 


0.7(5) 


2.3(7) 



TABLE III: Fits to (1 - qi)c with q^in = and q^ax = 0.4 
for periodic boundary conditions. The three groups of data 
refer, from top to bottom, to the three fitting functions I, II, 
and III in Eq. (1231 . respectively, and Form IV in Eq. 12511 
with d-ds = 0.43. 
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value agrees with the estimate /i = 0.44±0.02 of PY from 
the scahng of the spin-overlap but incompatible, within 
the statistical error bars, with the result /i = 0.63 ± 0.03 
for free be. We will return in Section IIII Di on the origin 
of the discrepancy between free and periodic be. The 
inset of Fig. El shows that, with these values of /i and 
d — ds, the scaling form for {1 — qi)c, Eq. H14|l . is also 
well satisfied. Finally, we verified that, if c? — ds = 0.43, 
the unrestricted average (1 — qi) satisfies scaling, giving 
/i = 0.45 ± 0.02 in agreement with the estimate from B. 

Combining Eqs. (|40|l and H41|l . we obtain the estimate 
of 9' for periodic be: 

e' = n- [d- ds) = -OM±Om (periodic be) . (42) 

This is compatible with zero and, within the error bars, 
incompatible with the value 9' = 9 ~ 0.2, where 9 char- 
acterizes the energy of domain walls induced by bound- 
ary condition changes. A scenario in which 9' = and 
d~ds > is consistent with the TNT picture. Finally, we 
note that, although our analysis of the PY data uses dif- 
ferent quantities to extract exponents, our results agree 
with those given by PY. 

D. Discussion and summary of the results 

In the previous sections, we have seen that for both free 
and periodic be, the analysis of all the different observ- 
ables considered gives consistent results for the exponents 
d — ds and 9' under the assumption of minimal correc- 
tions to scaling. However, while the results for d — ds for 
free and periodic be agree with each other, the results 
for 9' apparently do not, having found 9' — —0.01 ± 0.03 
for periodic be and 9' = 0.19 ± 0.06 for free be. Since 
9', like d — ds, should not depend on the type of bound- 
ary conditions, the discrepancy must be due to different 
corrections to scaling for the two boundary conditions. 

Therefore, it is important to analyze further the correc- 
tions to scaling. First, we recall that the scaling plots for 
the quantity B in Figs. (free be) and 1101 fperiodic be), 
from which we have determined /i (and hence 9'), are 
obtained by imposing that the whole data set (namely 
all values of e and L) satisfies scaling with corrections to 
scaling smaller than the statistical errors, which are less 
than 1%. Clearly, this is a very stringent requirement. If 
we relax this requirement, allowing some corrections to 
(simple) scaling, we can accommodate a larger range of 
values for fj.. 

This is shown in Fig. 1111 which gives a scaling plot for 
free be, analogous to Fig. |5lbut assuming the value /i — 
0.42 determined from periodic be. The polynomial fitting 
curve was obtained by excluding from the fit the data 
points for L — A and 6. One can see that for larger L the 
data collapse reasonably well on the curve. The deviation 
of the L = 4, 6 data from the curve, less than 10%, is a 
measure of the corrections to scaling. Therefore, we see 
that corrections to scaling of less than 10% for the two 
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smallest sizes are sufficient to remove the discrepancy in 
the value of 9' between free and periodic be. We verified 
that the also the other quantities considered, namely (1 — 
g/)c, i?, {I — qi), (1 — g), can be fitted in a similar way. 

We also tried the converse operation, namely a scaling 
plot of the data for periodic be but using the value fi — 
0.63. We find that one can get a relatively good data 
collapse excluding the sizes L = 3, 4, and 5, which deviate 
from the scaling curve by less than 10%. However, now 
the data for a given e approach the scaling curve from 
the right side instead of from the left side as in Fig. 1111 
but since they have an upward curvature, the correction 
to scaling should change curvature twice as L increases, 
which is not very plausible. Hence, we believe that it is 
more natural to conclude that the correct value of ^ is 
closer to 0.42 than to 0.63, namely that corrections to 
scaling are smaller for periodic than for free. 

Indeed, in general, it is reasonable to expect that cor- 
rections are larger for free be, because these be have a free 
surface on which lie a fraction of sites which is quite sub- 
stantial for moderate sizes. In Fig. ^] we plot together 
the (1 — qi)c data of Fig. 0| and El for free and periodic 
be. The data for free be lie significantly below those for 
periodic be, indicating that the surface of the excitations 
is smaller for free be. For periodic be, the domain wall 
has to "bend" to return to the same point on the "top 
surface" as it had on the "bottom surface" . This may be 
the source of the extra surface area. 
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Under the hypothesis d — dg ~ 0.44, Fig.^]then shows 
that the corrections to asymptotic scahng are larger for 
free be, since the free be data show a more marked devi- 
ation from the asymptotic e-independent behavior, and 
display a larger curvature. This is further indication that 
free be have larger corrections. 

Evidence that free be have larger corrections was also 
found recently in Monte Carlo simulations^^, where some 
evidence was observed that the free be data might have 
a crossover from droplet-like to either TNT- or RSB-like 
behavior at large sizes. 

Incidentally, note that if RSB is the correct asymptotic 
picture and the L —^ oo limit of (1 — 9;)c is the same for 
periodic and free be, then Fig. ^1 would indicate that 
the corrections are smaller for free be (since the data 
are closer to their non-zero asymptotic value) which is 
not very plausible. Note, however, that we do not have 
an argument why in the thermodynamic limit (1 — qi)c 
should be independent of boundary conditions. 

To summarize the first part of the article, we have an- 
alyzed several quantities for periodic and free be. For 
both types of boundary conditions, all the data are well 
described by a general scaling picture involving only two 
scaling exponents, d — dg and 9\ with only small correc- 
tions to scaling. Some observables show significant cor- 
rections to asymptotic scaling, which are larger for free 
boundary conditions. Fitting this scaling picture to our 
data, we obtain comparable values of d — ds for periodic 
(0.43 ±0.02) and free boundary conditions (0.44 ±0.03). 

By imposing that corrections to scaling are less than 
the statistical errors of 1%, for periodic boundary condi- 
tions we obtain 6' ~ 0, which fits well the TNT scenario 
{d — ds > 0, 6' = 0), while for free boundary conditions 
we obtain 6' = 0.19 ± 0.06, which fits well the droplet 
picture {d — ds > 0, 9' > 0). By relaxing this require- 
ment and allowing larger corrections to scaling of order 
10%, the data for free be can be also fitted by a sce- 
nario with 6*' ~ 0. Therefore the data for free be are 
also consistent with the TNT picture provided moderate 
corrections to scaling are allowed, larger than those for 
periodic be. We have also provided direct evidence that 
indeed free be have larger corrections to scaling. 

Data for the box overlap for free boundary conditions 
indicates smaller corrections to asymptotic scaling, which 
is reasonable since the box is away from the surface, and 
are consistent with the scenario described above. 

For both free and periodic be, the data are also fitted 
well by the RSB picture (d - 4 = 0, 9' = 0), but only 
if we allow very large corrections to scaling. In this case, 
the good scaling behavior we observed for all the observ- 
able considered would only be a finite size artifact, and 
would disappear at larger sizes. To test this possibility, 
large system sizes will be needed. 

This concludes the first part of the paper, dedicated to 
the physical results. In the second part, we will describe 
the branch and cut algorithm employed, and analyze its 
performance in our computations. 
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IV. BRANCH AND CUT ALGORITHM 

Branch and cut is, to our knowledge, the fastest exact 
method for determining ground states of spin glasses in 
three dimensions. To apply this technique, we transform 
the problem of minimizing the Hamiltonian in Eq. into 
a standard combinatorial optimization problem known as 
the maximum cut problem. (For a detailed description of 
optimization and related topics, see Ref. |^.) Consider 
the interaction graph G = {V^ E) associated with the spin 
glass Hamiltonian, where G contains vertices 1 , . . . , L"^ e 
V associated with the spin sites and edges iij) E E with 
weight Cij — —Jij associated with the couplings. 

Given a partition of V into two sets, W C V and its 
complement V \ W, the cut S{W) associated with W is 
defined as the set of edges with one endpoint, i say, in 
W and the other endpoint, j say, in 1/ \ W. In formulas, 
S{W) = {{ij) e E \ i e W,j e V\W}. The weight of 
a cut S{W) is defined as the sum of the weights of the 
cut edges X](ij)Gi5(VK) ^ m,aximum cut is a cut with 
maximum weight among all partitions W. It is easy to 
show that minimizing the Hamiltonian Eq.l^ is equiva- 
lent to finding a maximum cut in G, see Ref. Il9t If we 
know a maximum cut with node partition W and V\ W, 
the corresponding ground state spin configuration can be 
read off by assigning the value up to the spins in W and 
down to the spins in 1/ \ W, or vice versa. 

The branch-and-cut algorithm solves the maximum cut 
problem through simultaneous lower and upper bound 
computations. By definition, the weight of any cut gives a 
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lower bound on the optimal cut value. Thus, we can start 
from any cut and iteratively improve the lower bound us- 
ing deterministic heuristic rules (local search and other 
specialized heuristics, see RefiSli for details). How do 
we decide when a cut is optimal? This can be done by 
additionally maintaining upper bounds on the value of 
the maximum cut. Upon iteration of the algorithm, pro- 
gressively tighter bounds are found, until optimality is 
reached. 

Since the availability of upper bounds marks the differ- 
ence between a heuristic and an exact solution, we now 
summarize how the upper bound is computed (for more 
details, see Ref.lsil) To each edge (ij) we associate a real 
variable Xij and to each cut S{W) an incidence vector 
T^s{w) g components xfj'^'' associated to each 

edge (ij), where xtj^^ = 1 if (ij) G 6(W) and xtj^^ = 
otherwise. Denoting by Pc{G) the convex hull of the in- 
cidence vectors, it can be shown that a basic optimum 
solution—' of the linear program 

max{ ^ J.jX^j I X G Pc{G)}. (43) 

(ij)£E 

is a maximum cut. In order to solve H43|) with linear 
programming techniques we would have to express Pc (G) 
in the form 

Pc{G) = {x eM.^ \ Ax <b,0 < X < 1} (44) 

for some matrix A and some vector b. Whereas the exis- 
tence of A and b are theoretically guaranteed, even sub- 
sets of Ax < b known in the literature contain a huge 
number of inequalities that render a direct solution of 
(|43|l impractical. 

Instead, the branch-and-cut algorithm proceeds by op- 
timizing over a superset P containing Pc (G) , and by it- 
eratively tightening P, generating in this way progres- 
sively better upper bounds. The supersets P are gener- 
ated by a cutting plane approach. Starting with some P, 
we solve the linear program niax{^^j^-^g^ JijXij \ x G P} 
by Dantzig's simplex algorithmS^. Optimality is proven 
if either of two conditions is satisfied: (i) the optimal 
value equals the lower bound; (ii) the solution vector x 
is the incidence vector of a cut. 

If neither is satisfied, we have to tighten P by solving 
the separation problem. This consists in identifying in- 
equalities that are valid for all points in Pc{G), yet are 
violated by x, or reporting that no such inequality ex- 
ists. The inequalities found in this way are added to the 
linear programming formulation, obtaining a new tighter 
partial system P' d P which does not contain x. The 
procedure is then repeated on P' and so on. 

At some point, it may happen that (i) and (ii) are not 
satisfied, yet the separation routines do not find any new 
cutting plane. In this case, we branch on some fractional 
edge variable Xij (i.e. a variable Xij ^ {0,1}), creating 
two subproblems in which Xij is set to and 1, respec- 
tively. We then we apply the cutting plane algorithm 
recursively for both subproblems. 



V. PERFORMANCE OF THE BRANCH AND 
CUT ALGORITHM 

In this section we study the performance of our cur- 
rent implementation of the branch and cut algorithm, in 
particular the dependence of the number of computer op- 
erations on system size. The results for size L = 12 were 
obtained with a more efficient version of the code, so per- 
formance for this size cannot be compared with that for 
the smaller sizes. Hence, in this section, we shall just 
consider sizes up to L = 10. 

Finding the ground state of the Hamiltonian Eq. 
in three dimensions is an A/'P-hard problem.^^, and all 
known algorithms to solve this class of problems require 
a number of operations that grows exponentially on the 
size of the input, in the worst case. 

However, depending on the problem, the number of 
operations for typical instances (for the spin glass prob- 
lem, an instance is a realization of the random couplings, 
or sample), can grow considerably more slowly than the 
worst case exponential behavior. Furthermore, the num- 
ber of operations can vary significantly from one instance 
to another. It is therefore useful to investigate experi- 
mentally the performance of the algorithm for typical in- 
stances, in order to try to extrapolate the computational 
resources necessary to go to larger sizes, and possibly to 
identify which parameters of the problem affect most the 
performance. De Simone et al^ measured the average 
CPU time used by the branch and cut algorithm to find 
the ground state of the two-dimensional ± J spin glass 
with periodic be, up to L = 70, showing that the average 
CPU time was approximated by a function proportional 
to L^. 

Here we analyze the performance of the branch and 
cut algorithm for the three-dimensional spin glass with 
free be and Gaussian couplings. In order to do this, we 
first need a good measure of the performance. For a 
complex algorithm such as branch and cut, a simple and 
absolute measure of the number of operations is not avail- 
able. Two possible measures are the CPU time and the 
number of linear programs solved during the run of the 
algorithm. In Table Hvl we summarize the average run- 
ning time needed for calculating an unperturbed ground 
state for the different system sizes. 

The CPU time is not an accurate measure since it de- 
pends on the machine architecture and load. Further- 
more, our computations were carried out on several dif- 
ferent machines, so the CPU time is not useful here. We 
take instead the number of linear programs solved, Up, 
because (i) it is a well-defined and machine independent 
quantity; (ii) we have observed that about 95% of the 
time is spent in solving linear programs; (iii) for a fixed 
system size, Up correlates strongly, and almost linearly, 
with the CPU time. This is shown in Fig. 1131 which plots 
the CPU time versus Up for 1000 randomly generated 
samples with L = 10, computed on the same machine. 
Note that since the size of the linear programs is also 
growing with the system size, the CPU time per linear 
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FIG. 13: Scatter plot of the CPU time to find the unper- 
turbed ground state (e = 0) versus the corresponding number 
of Hnear programs solved {rip). Each point represents a ran- 
domly generated sample with L — 10. All the computations 
for this set of samples were run on the same machine. The 
dashed line indicates a linear behavior. 



L 


mean CPU time per sample 


4 


0.065 


6 


0.662 


8 


10.11 


10 


338 



TABLE IV: Mean CPU time per sample in seconds for the 
calculation of the unperturbed ground state, averaged also 
over the different machines. 



program increases strongly with L: the average (resp. 
median) CPU time goes from 0.00770 (resp. 0.833) sec- 
onds for L = 4 to 0.833 (resp. 0.784) seconds for L = 10. 

Hence, Up severely underestimates the rate at which 
the number of operations increases with L. 

From Fig. ^1 we also note that the distribution of 
Up (and CPU times) is very broad, extending over three 
orders of magnitude. The histogram distribution of Up 
for different system sizes is shown in Fig. 1141 In addi- 
tion to shifting to larger Hp, the distribution broadens as 
L increases. Also, there is some evidence of a double- 
peak structure. For L — 10, we verified that the peak at 
smaller Hp corresponds to samples that could be solved 
without branching, while the other peak corresponds to 
samples where branching was necessary. Since in each 
branching step the number of subproblems to be solved 



FIG. 14: Histogram of the number of linear programs solved 
by the branch and cut algorithm to find the unperturbed 
ground state for different system sizes. 



doubles, the number of linear programs increases rapidly 
and the second peak is at significantly larger Up. 

In order to identify which parameters of the problem, 
in addition to the size, affect the performance, we ask 
whether Up correlates with the physical observables we 
measure. No significant correlation was observed with 
the ground state energy. Fig. El plots^ (logio ^p) for 
the unperturbed ground state (e = 0) and L = 10 versus 
the overlap between this state and the perturbed ground 
state with e/r = 4. We observe a distinct correlation 
between Up and q: for small q, more linear programs 
are needed than for large q. The figure shows that the 
typical number of linear programs is close to an order of 
magnitude larger if (7 ~ than if g ~ 1. We observed a 
similar correlation for other values of e as well, and also 
between the CPU time and q. Again, the distribution of 
Up is quite broad as shown by the data for the standard 
deviation of logj^g Up in Fig. [T51 

In order to quantify how the correlation between Up 
and q changes with the system size, we show in Fig. II6I 
the average and median of Tip clS cL function of Nb, as well 
as the conditional averages of Up restricted to samples 
with large {\q\ > 0.9) and small {\q\ < 0.1) overlap. We 
take the number of bonds, Nb^ as a measure of the input 
size, since the maximum cut problem is specified in terms 
of the edge variables in the graph. From Fig.^jwe see 
that, first, all measures show an approximately exponen- 
tial increase with Nt, with corrections for small Nf,, and 
second, the difference between the conditional averages 
with small and large q seems to increase with the system 
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FIG. 15: The circles are a plot of (log^Q np), where is the FIG. 16: Average rip, median rip, and conditional averages 
number of linear programs solved to compute the unperturbed 
ground state 5''', versus the overlap between S'"' and the 
perturbed ground state S*-"'. The data is for e/r = 4 and 
the samples were selected from a set of randomly generated 
samples with L = 10, in such a way that the same number 
of samples is plotted for each consecutive q interval of length 
0.1, in order to sample equally all regions of g. The triangles 
show the standard deviation, among samples, of logj^Q Tip as a 
function of q. 



of rip restricted to \q\ < 0.1 and to |g| > 0.9, as a function 
of the number of bonds A^i,. The data for rip are for L = 10 
and e = (unperturbed ground state), and q is the overlap 
between the e = and e/r — 4 ground states. 



size, and is about one order of magnitude for L = 10. 

A qualitative difference between samples with small 
and large overlap is that samples with a small \q\ have a 
rougher "energy landscape" , namely states with an en- 
ergy close to the ground state energy yet a spin configura- 
tion very different from the ground state. It is then intu- 
itively clear why one would observe a correlation between 
q and the running time for a stochastic algorithm employ- 
ing local search heuristics, such as simulated annealing, 
since when the algorithm encounters one of these config- 
urations with small overlap, it must retrace its steps by 
a large amount. 

For the branch and cut algorithm, the reason for the 
correlation between Up and q is less obvious, but some 
insight is provided by an analysis of "reduced cost fix- 
ing". This is a feature of the branch and cut algorithm 
speeding up the computations. In every iteration of the 
algorithm, reduced cost fixing gives us a sufhcient condi- 
tion to decide which variables (corresponding to the edges 
in the graph) have already attained their optimal value. 
Therefore, we can fix the variables with "optimal" status 
to their current value for all the subsequent iterations of 
the algorithm, resulting in less overall computational ef- 
fort. The more variables that can be fixed, the faster the 



algorithm is in practice. 

Since the samples with small overlap have "almost op- 
timal" solutions with spin configurations very different 
from the ground state, a smaller number of variables can 
be fixed. Here we do not have the "correct" edge values 
available until the end. As an example, we checked that 
for L = 10 and e = t, for 100 randomly chosen samples 
with small overlap {\q\ < 0.1), in average 409 ± 39 of the 
2700 edge variables could be fixed in the first sub prob- 
lem, i.e. before branching takes place. In contrast, for 
100 randomly chosen samples with big overlap {\q\ > 0.9), 
921 ±34 of the edge variables could be fixed in the first sub 
problem, about twice as many. Of course, the less vari- 
ables that can be fixed in the first sub problem, the more 
overall branching is necessary, resulting in more overall 
computational effort for samples with small overlap. 

A consequence of the the broad distribution of the 
CPU time and of its correlation with the physical ob- 
servables of interest, is that a cutoff in the CPU time 
produces a systematic error in these quantities. One has 
therefore to ensure that the cutoff is large enough so that 
the systematic error is smaller than the statistical error. 

It is interesting to try to extrapolate the running time 
needed to deal with larger sizes. The average CPU time 
in Table IWI varies approximately as ~ exp(aA^b) with a 
somewhere between 0.0024 and 0.003. Extrapolating to 
L = 14 {Nf, = 7644), this gives an average CPU time 
of around lO'^*^ seconds per sample, which is clearly 
very demanding. Furthermore, memory limitations will 
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set in before we can reach this size. Again, note that 
Up increases much more slowly with N},. The data for 
l^l < 0.1 in Fig llfil for example, vary approximately as 
~ exp(Q:iVt,) with a smaller a around 0.0017, showing 
that the dominant limiting factor is the solution of the 
linear programs. Note that the program used for L = 12 
is significantly faster than that used in this extrapolation. 
This long extrapolated running time gives us further mo- 
tivation to continue our research on the improvement of 
this algorithm. 



d — ds and 9' , and how the values of these parameters 
predicted by the various physical pictures proposed for 
the spin glass phase fit our data. Our conclusions have 
been summarized at the end of Section IlII Dl 

In the second part of the paper, we have analyzed the 
performance of the branch and cut algorithm, finding 
that the performance is worse when there is a low energy 
excited state close in energy to the ground state but far 
away in configuration space, and have given a quantita- 
tive analysis of this effect. 



VI. CONCLUSIONS 

Using an exact "branch and cut" optimization algo- 
rithm, we have studied the large-scale, low-energy exci- 
tations in the Ising spin glass in three dimensions with 
free boundary conditions, and compared the results with 
those obtained earlier by PY for periodic boundary con- 
ditions. 

In the first part of the paper, we have discussed in de- 
tail how the whole set of observables analyzed is fitted by 
a general scaling picture characterized by two exponents. 
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